The Time-Dependent Local Density 
Approximation for Collective Excitations 

of Atomic Clusters 

G.F. Bertsch 

Physics Department and Institute for Nuclear Theory 
University of Washington, Seattle, WA 98195 USA 
E-mail: bertsch@phys.washington.edu 
and 

K. Yabana 

Department of Physics, Niigata University 
Niigata, Japan 
E-mail:yabana@ carrot, sc. niigata-u.ac.jp 



Abstract. We discuss the calculation of collective excitations in atomic clusters using 
the time-dependent local density approximation. In principle there are many formula- 
tions of the TDLDA, but we have found that a particularly efficient method for large 
clusters is to use a coordinate space mesh and the algorithms for the operators and the 
evolution equations that had been developed for the nuclear time-dependent Hartree- 
Fock theory. The TDLDA works remarkably well to describe the strong excitations 
in alkali metal clusters and in carbon clusters. We show as an example the benzene 
molecule, which has two strong features in its spectrum. The systematics of the linear 
carbon chains is well reproduced, and may be understood in rather simple terms. 

I INTRODUCTION 

The time-dependent local density approximation (TDLDA) is a powerful tool to 
calculate the quantum mechanical motion of electrons in condensed matter systems. 
In my talk I want to describe to you the numerical techniques we use [1], which 
were borrowed from nuclear physics [2] . I then want to show you a survey of some 
of the results, concluding with new work on the behavior of electrons in elongated 
chains. The time-independent LDA is known as Density Functional Theory and 
is well-established in condensed matter physics [3] as a predictive ab initio theory 
that is practical for structures beyond the range of quantum chemistry methods 
[4-7]. The equations of the LDA are easy to write down. The energy function is 
given by the expression 
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where the fa are the single-particle wave functions and p(r) = ^Zi\<fii{r)\ 2 is the 
density. The potentials included in the energy function are V ion , the ionic poten- 
tial, e 2 /\r — r'\, the direct Coulomb interaction between electrons, and v ex (p), the 
local-density approximation to the exchange-correlation energy. The static LDA 
theory is obtained by minimizing this energy function, requiring only that the fa 
be orthonormal. That gives the Kohn-Sham equations for the wave functions fa. 
The TDLDA equations are very similar, with the "energy" (actually the Lagrange 
multiplier) in the Kohn-Sham equation replaced by the time derivative ihd t , 

thd t fa(r,t) = Hfa(r,t). (2) 

with 
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This is almost all I want to say in general. One always uses pseudopotentials for 
the ions [8,9], to eliminate the core electrons with their large energy scales. The 
LDA has well-known difficulties in describing the electronic excitations of systems, 
attributable to the simplified treatment of exchange. However, the collective motion 
is rather insensitive to the nonlocality of the exchange, and the TDLDA is much 
more reliable than one would expect from looking at energy gaps. 



II NUMERICAL 

For most applications it is sufficient to consider small deviations from the static 
solution, and then there are a number of approaches to solve the small-amplitude 
TDLDA equations. Nuclear physicists are most familiar with the configuration rep- 
resentation, which I'll just mention briefly. Here one represents the time-dependent 
wave function in the basis of the static solutions, 
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The equation of motion is Fourier transformed to frequency, and the matrix equa- 
tions satisfied by the particle-hole amplitudes dij(u) are the RPA equations of 
motion. This method becomes inefficient for large systems, because the the di- 
mensionality of the matrix grows as the product of the number of particle and 
hole orbitals. In the more efficient methods, the dimensionality increases linearly 
with the number of particles. One method that is very popular is to calculate the 



density-density response function. This is the density change induced at a point r 
by a density perturbation at a point r' . Because of the local density approximation, 
one can calculate the full response from the response of non-interacting electrons. 
The number of mesh points to represent the response increases as the volume of 
the system, and thus only linearly with the number of particles, making it more 
efficient for large systems. The method has been used extensively for systems with 
spherical symmetry [10,11]; recently it has been applied also to clusters in three 
dimensions [12]. 

We use another method whose dimensionality also scales linearly with the system 
size. This is the straightforward solution of the time-dependent equations of motion. 
The number of variables is the number of points at which wave function amplitudes 
are represented. This can be of the order of 50,000 for some of the examples we 
consider, so some care must be taken to use efficient numerical algorithms. We 
have simply adopted the technique of ref. [2] who investigated the TDHF theory of 
nuclear dynamics. 

We use a uniform mesh in coordinate space to represent the wave function. The 
shape of the gridded volume can be spherical, cylindrical, or rectangular, depending 
on the cluster under consideration. If the equations were linear in 0, they could 
efficiently integrated by using the Taylor series expansion of the evolution operator, 
4>(t) = exp(—iHt)(f) = 1 +tH(f) + t 2 H 2 <f)/2 + ... In practice we can ignore the change 
of H dependence of H on for some small time interval r, and use the Taylor series 
integrator over that interval. We typically use fourth order: 
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Considerable care must be given to defining the H to propagate in the above 
equation. An important consideration is that the algorithm conserve energy to a 
very high accuracy. This can be assured with an implicit definition of H in terms 
of the densities pit) and pit + r) at the beginning and end of the interval r. The 
definition is 

Yy J L 2m p(r,t + r) - p(r,t) J w w 

In practice we do not calculate this exactly but rather estimate the Hamiltonian at 
t + r/2 using a predictor-corrector method. In our algorithm, first the potential is 
extrapolated forward by half a time step with a 3-point difference formula to get a 
trial potential V tr i a i{t + r/2). The time integration is then performed using V tr i a i 
to get a trial density ptriai{t + t), and the density at time t + r/2 is estimated as 
pit + r/2) = {pit) + ptriai)/2. The potential obtained with this density is used to 
once more integrate the wave function from t to t + r. 

To generate the response of the system, we perturb the ground state wave func- 
tion and measure the response as a function of time. The initial wave functions are 
taken as <f>i{f, 0) = e lkz (f)°{r) in which (f)° represent the static Kohn-Sham solutions. 



We are interested in the optical response, requiring the perturbing parameter k to 
is small. This ensures also that the response will be in the linear regime. The real 
time evolution of the dipole moment is obtained as z(t) = J2i < 4 > i{t)\ z \4 > i{t) >> 
and its Fourier transform in time gives the dipole strength function. 

Summarizing, our algorithm to integrate the TDLDA equations has as numerical 
parameters the spatial mesh size, Ax, the number of mesh points M, the time step 
r, and the total length of time integration T. The actual values for the different 
systems will be given below. 

Ill PHYSICS 

Our interest is in the electronic excitations, and we will leave aside completely 
the problem of determining the ionic geometries of the clusters. In classical electro- 
magnetic theory, the dipole strength in spherical metal clusters is concentrated in 
the Mie resonance. The resonance may be sharp or not depending on the cluster, 
and it may be present also in clusters of different shapes. 

A Alkali metals 

The simple, sharp Mie resonance behavior is reproduced by jellium model of 
alkali metal clusters [10], which however neglect the perturbing effects of the ionic 
cores. It is therefore interesting to see to what extent the ionic cores affect the 
physics. I will first consider a system where the jellium approximation is remarkably 
good. This is the for clusters of sodium atoms. An example of a calculation is the 
sodium cluster Na^. This was chosen because 147 atoms can be placed in a 
icosahedral structure, which is as close a one can get to spherical structure with 
discrete atoms. The numerical parameters of the calculation are: Ax = 0.8 A, 
r/h = 0.01 eV" 1 , M = 28000 and N t = 5000. Because the valence electron is in an 
s-wave, and the short-distance part of the wave function has been smoothed away 
by the pseudopotential, a rather coarse Ax is sufficient. The total time integration 
is rN t /h = 50 eV -1 which in principle allows structures on the scale of 0.02 eV 
to be resolved. Fig. 1 shows the response, compared with the spherical jellium 
model. Both descriptions show a sharp peak at just over 3 eV excitation, and some 
background extending in both directions. The two theories are quite close. This 
shows that the pseudopotential in sodium has little effect. Experimentally, sodium 
clusters show a fairly sharp resonance, but the energy is lower by about 15%. 

The next case, lithium, shows a more dramatic difference between the jellium 
model and the full pseudopotential TDLDA. Fig. 2 shows the response of Lijj 7 , 
comparing the full pseudopotential LDA with the jellium model. The Mie resonance 
has been shifted down by the pseudopotential interactions and broadened. The 
reason for the difference compared to sodium is not hard to find [14]. The scattering 
effect of the pseudopotential is much larger in the case of Li, and the scattering 
gives the electrons an effective mass larger than one. In the jellium model the ion 
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FIGURE 1. Dipole response of a sodium cluster, Na+ 47 , comparing the spherical jellium model 
with the full pscudopotential TDLDA calculation. 
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FIGURE 2. Dipole response of lithium clusters, comparing the jellium model, the full TDLDA 
calculation, and experiment [13]. 



behaves as a uniform positive charge density filling space within the lattice. The 
scattering effect of the pseudopotential may be seen by examining the scattering 
from the difference between the pseudopotential and a uniform positive charge. 



Thus, we are led to examine the scattering from the difference potential, 

AV-W + r<r, 
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When we do this we find that both the s- and p-wave scattering phase shifts for 
sodium are small, showing that effects of the pseudopotential are weak. On the 
other hand, in the case of Li, the s-potential is repulsive but the p-wave potential is 
attractive. This implies that there will be a large backward angle scattering, making 
the single-particle wave functions rather different than the jellium wave functions. 
We believe this aspect is responsible for the broadening of the Mie resonance in Li. 
Another consequence of the two potentials is that the effective forward scattering 
potential will be momentum-dependent, becoming more attractive as the momen- 
tum increases and the amount of p-wave increases. This lowers the Mie resonance 
frequency. 

B Carbon structures 

Carbon is more difficult because the p- valence electrons are rather tightly bound 
and require a finer mesh. In practice, a mesh size of Ax = 0.3 A is adequate to 
calculate the electronic response of carbon structures, given their geometry. As a 
first test, we examine the response of the benzene molecule, which has been well- 
studied up to 30 eV with synchrotron photon sources [15,16]. For the calculation, 
we used a pseudopotential for the hydrogen atoms as well as the carbon atoms, 
and used the actual geometry of the molecule. The grid of mesh mesh points 
was limited to points in a sphere of radius 7 A. This requires about 50,000 mesh 
points. Finally, because of the tighter binding of the p-waves, a time step a factor 
of 10 smaller than that for the alkali metals was required, r = 0.001 eV -1 . The 
resulting response is plotted in Fig. 3, showing the Fourier transform of the real- 
time response computed over a time interval of 30 h/eV. This is compared to the 
strength function measured in ref. [15], which an adjusted normalization (which 
was arbitrary in ref. [15]. The main features of the response are a sharp peak 
at 6.9 eV, and a broad bump centered just below 20 eV. The sharp peak is very 
well reproduced by theory. It is associated with a collective resonance of the n 
electrons. The calculated strength of the resonance / is =1.1 in the usual sum rule 
units (i.e. the effective number of electrons). The experimental strength can be 
extracted from the data of ref. [16], but the peak shape has wings that are not 
simple Lorentzian tails. Including the right-hand wing up to 8.3 eV, but not the 
left-hand shoulder, the measured strength is / = 0.9, which we consider rather 
good agreement. The large broad bump is associated with the a electrons. We 
reproduce its gross structure. But note that the theoretical distribution has much 
finer structure than is observed experimentally. 
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FIGURE 3. In-plane dipole response of benzene molecule, comparing TDLDA with experiment 
[15]. 



Our next example is Ceo- Our TDLDA calculation is shown in Fig. 4, together 
with another TDLDA calculation using the matrix RPA method [18]. In our cal- 
culation, the 7r electron excitation appears as a single peak at 7 eV, carrying a 
strength of / = 9. Experimentally, the 7r electron transitions are split into several 
peaks, with a total strength up to 7 eV of / = 6. There are no absolute measure- 
ments of the strength above 7 eV, but photoionization measurements exist [17] that 
show a broad a resonance, as in benzene. The (unnormalized) data is shown in the 
figure with diamond symbols. The polarizability of has some interest, because 
the simplified tight-binding model gave a number much smaller than experiment, 
a 45 A 3 compared to the experimental value a 85 A 3 . The polarizability may 
be extracted from the dipole strength function by the integral 

_ e 2 r df dE 
01 ~ mJ dE E 2 

The result is 80 A 3 , in satisfactory agreement with experiment. 

C Carbon chains 

The important role of geometry in the dipole response is beautifully illustrated 
by the collective ir transition in carbon chains and rings. For these calculations 
[19], we guessed at the geometry, fixing the nearest neighbor distance at 1.28 A for 
all cases. This is the average LDA equilibrium distance for large rings or chains. 
The atoms are in a line on the chains and form a circle on the rings. As in the 
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FIGURE 4. Response of C 60 in TDLDA. The dashed curve is another TDLDA calculation, 
using the particle- hole matrix formulation [18]. 



other examples, there is a strong collective 7r transition along the axis of the chain 
or in the plane of the ring. The energies of the excitations are plotted in Fig. 5 
as a function of the number of atoms. There is experimental data showing the 
existence of excitations at the predicted energies for chains [20]. Unfortunately, 
the strengths of the transitions were not measured, so it is not known whether the 
observed transitions really correspond to the calculated collective excitations. The 
stable form of carbon for cluster numbers greater than ten or so is believed to be 
the ring configuration, but no data is available on this form. 

Notice that the excitation energies have a strong but smooth dependence on the 
chain length. The variation is entirely explicable in terms of simple physics, not 
requiring at all the detailed TDLDA calculations. This is the theme of the next 
section of my talk. 



IV COLLECTIVE ENERGETICS IN ATOMIC CHAINS 

The behavior of electromagnetic resonances on infinitely long wires is known from 
classical electromagnetic theory. The dispersion formula for the one-dimensional 
plasmon on a long wire reduces to the following expression in the long-wave length, 
thin wire limit [23,24]. 

2 4vrn e e 2 2 1 

u = q log— (7) 

m qa 



where q is the reduced wave number of the plasmon, n e is the density of electrons 
per unit length, and a is the radius of the wire. For a finite wire, the lowest mode 
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FIGURE 5. Excitation energies of TDLDA collective states and single-particle excitations in 
carbon chain and rings. The crosses mark experimentally observed excitations [20-22]. 



would have a q varying inversely with the length of the wire L. Thus the lowest 
mode would behave as 



uj ~ C (8) 

Lj 

This behavior can be extracted from a more quantum approach is the polariz- 
ability estimate of the collective frequency 



2 h 2 e 2 N e 

ui = (9) 

ma 

where N e is the number of active electrons and a is the polarizability. This formula 
is derived from the ratio of sum rules, and N e may be identified with the oscillator 
strength / associated with the transition. For the linear carbon chains, the transi- 
tion is associated with n electrons and and the number of them in the chain C n (n 
odd) is 

N e = 2n- 2. 

The TDLDA calculations confirm that this is satisfied at the level of 20% accuracy. 
The polarizability is harder to estimate. In ref. [19] we model the polarizability 
assuming that the chain acts as a perfectly conducting wire, of some fixed trans- 
verse dimension a and a long length L. Then one can show that for large L the 
polarizability is given by 

L 3 
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FIGURE 6. Comparison of analytic formulas with the energy systematics of carbon chains. 
Solid line, eq. (10); dashed line, eq. (11). 



The dependence on the number of atoms in the chain follows form L ~ n. Then 
inserting the above in eq. (9), we find the large-n behavior 

y/\ogn 

n 

This is plotted in Fig. 6 as the solid line, fitted to the TDLDA result at n — 14. 
We see that the general trend is reproduced, but the asymptotic behavior is only 
realized in long chains. 

Another classical view is to compare with the polarizability of a conducting 
ellipsoid, which can be calculated analytically [26]. The formula is 

where e is related to the ratio of short to long axes, R±/R 

This has been applied to the longitudinal mode in fullerines [6]. Its asymptotic 
behavior is given by eq. (10) above. Remarkably, this formula gives an excellent 
fit, treating the two parameters as adjustable. This is shown as the dashed curve 
in Fig. 5. 

We finally discuss the relative frequencies of the modes in chains and rings. One 
expects that for C n in the form of rings, the number of n electrons is 2n, roughly 
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the same as chains. However, the polarizability will be quite a bit smaller because 
of the more compact geometry. The ratio of polarizabilities is about a factor of 
two, giving a prediction of a 40% higher collective resonance. In fact the TDLDA 
calculations shows that the sum rule / is also higher for ring configurations, and 
the actual resonance is at about twice the energy of the chain. This may be seen 
in Fig. 5. 

V WHAT NEXT? 

We have discussed carbon, but from a computational point of view the next 
group IV element, silicon, would be very similar. The next challenge for numerical 
TDLDA would be the group IB metals, the so-called coinage metals. These have 
a single valence electron in an s shell, like the alkali metals, but there is a nearby 
closed d shell that cannot be neglected in the response. Beyond that, there are all 
rf-shell metals which exhibit broad responses, hardly showing any trace of the Mie 
resonance. But this remains for the for the future. 
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